Interplay between the Ionic and Electronic Density Profiles in Liquid Metal Surfaces. 
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First principles molecular dynamics simulations have been performed for the liquid-vapor inter- 
faces of liquid Li, Mg, Al and Si. We analize the oscillatory ionic and valence electronic density 
profiles obtained, their wavelengths and the mechanisms behind their relative phase-shift. 
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X-ray reflectivity measurements on the surface of liq- 
uid metals and alloys, along with other techniques like 
diffuse scattering or grazing incidence diffraction, have 
shown the existence of layering in the ionic density 
profile i2iMi£i£iL£i Its origin is yet not clear and sev- 
eral reasons have been mooted. Rice and coworkers 910-11 
have pointed that the reason for surface layering in met- 
als is the interconnection between the ionic and electronic 
densities and that the abrupt decay of the electron den- 
sity at the surface induces an effective wall against which 
the ions, behaving like hard-spheres, stack. Recently, 
it has been suggested that surface layering is a rather 
universal phenomenon, although in most cases it is frus- 
trated by solidification^ 2 ^ therefore it only appears in 
systems whose melting temperature is very low compared 
with the critical temperature. 

The reflectivity experiments probe the total electronic 
density profile. Therefrom, the ionic density profile is 
derived by a superposition approximation where the the 
total electron density is taken as the sum of atomic-like 
electron densities around each of the nuclei in the sys- 
tem. Whereas this approach may be vindicated for the 
tightly bound core electrons, the case of valence electrons 
is more subtle. Early calculations of Lang and Kohn for 
semiinfinite step surfaces^ showed that the valence elec- 
tron density does not decay so abruptly, but displays 
some spill-out outside the surface; moreover inside the 
bulk its exhibits the so called Friedel oscillations. 

Computer simulations of liquid surfaces can evaluate 
directly the ionic density profile. But only a fully ab ini- 
tio method can deliver an electronic density selfconsistent 
with the discrete nature of the ions. Orbital based ab 
initio simulations are scarce (just for Sii£ and Nai&) due 
to the huge computational demands they pose. Orbital 
free ab initio simulations are less demanding, although 
still expensive, and recent orbital free ab initio molecu- 
lar dynamics (OFAIMD) calculations with 2000 and 3000 
particles have studied the surface properties of liquid Li, 
Na, Nao.3Ko.7 and Lio^Nao.e^*^ These studies showed 
that the superposition approximation produces a valence 
electronic density profile very similar to the fully selfcon- 
sistent one, except for the width of the interface due to 
the spill-out. 

This communication reports results for the liquid- 
vapor interface of liquid Li, Mg, Al and Si near their 
respective melting points. The calculations were per- 
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TABLE I: Thermodynamic states and simulation details. 



formed by the OFAIMD method where the forces acting 
on the nuclei are computed from electronic structure cal- 
culations, based on the density functional theory (DFT), 
which are performed as the ionic trajectories are gener- 
ated. In OFAIMD— the energy and forces acting on a 
system of N ions are computed from the ground state 
energy and density of the valence electrons, which inter- 
act with the ions through suitable local pseudopotentials. 
According to DFT the ground state electron density for a 
given ionic configuration minimizes an energy functional 
which is the sum of the electronic kinetic energy, the 
classical electrostatic energy, the energy of interaction 
with the ions through local pseudopotentials, and the 
exchange and correlation energy, for which we use the lo- 
cal density approximation. The keynote of the OFAIMD 
simulations is the use of an explicit, but approximate, ki- 
netic energy functional. Another basic magnitude is the 
local pseudopotential, v ps (r), describing the ion-electron 
interaction. It has been developed from first principles 
by fitting to the displaced electronic density induced by 
an ion immersed in a metallic medium. Further details 
on the method appear in references |20ll21| . 

The application of this formalism to bulk liquid Li, 
Mg and Al22i2I has provided an accurate description of 
their static and dynamic properties. Recent calculations 
for bulk liquid SiS^ have yielded a good description of 
the static structure factor (i.e. peak positions and am- 
plitudes as well as the shoulder at the high-g side of the 
main peak); moreover, the number of neighbors is around 
6, as in the experiment, while the diffusion coefficient is 
in the same range as other ab initio results. 

The simulations proceed as follows: given the ionic po- 
sitions at time t, the electron density is expanded in plane 
waves with energy less than a given cutoff, -E C ut, and 
the energy functional is minimized with respect to the 
plane wave coefficients, yielding the ground state elec- 



2 




-10 -5 



1 


A 


i 




AT V\ 





1.1 



1.2 
1 

0.8 



^ 3 



K 



*Rb 



Ba" 



Tl»Na 
Si Ga # »Sn 

It I 







Ba 



K 



Cs • 

Rb 



■5n Tl 

•fGa *Li 

It i I i oj I 



•Na 







1 



R WS (A) 



2 

r(A) 



FIG. 2: Wavelength of the oscillations in the density profiles 
as a function of the Wigner-Seitz radius and r 3 , the radius of 
a sphere which on average contains one electron. 



FIG. 1: Normalized ionic (thick line) and valence electronic 
(thin line) density profiles normal to the liquid-vapor inter- 
faces. The Mg, Al and Si data are shifted upwards by 0.3, 0.6 
and 0.9 units. The insets highlight the region near the surface 
to enhance the valence electronic density oscillations. 

tronic density and energy, and the forces on the ions, 
based on which the ionic positions and velocities are up- 
dated. For all the systems, the associated density pro- 
files were computed based on a sample of 20000 config- 
urations. The simulation setup consists of a periodically 
repeated supercell with dimensions L x L x L z , and a liq- 
uid slab of 2000 particles placed initially in the center of 
the supercell, occupying a volume consistent with the ex- 
perimental densities at the temperatures considered, and 
two free surfaces normal to the z axis. In all cases the 
distance, d, between the two surfaces is greater than 45 
A and the distance, 5 — L z — d, between the periodically 
repeated slabs is greater than 14 A, which are enough to 
guarantee the absence of unwanted interactions between 
the slabs and between the two surfaces of one slab. Table 
[3 summarizes these data and other simulation details. 

The longitudinal ionic density profiles were computed 
from a histogram of particle positions relative to the slab 
center of mass and the results are depicted in figure ^ 
All systems show stratification for around four layers into 
the bulk liquid, which agrees with the experimental ob- 
servations in other liquid metals. The wavelength of the 
ionic oscillations, A, shows a good scaling with the radii 
of the associated Wigner-Seitz spheres, Rws', however 
no clear relationship with electronic parameters, like the 
radii per electron, r s , has been found (see figure[21 where 
we also include data for other systems studied within the 
same method). This fact suggests that the ionic oscilla- 
tions are not induced by the Fricdcl oscillations in the 
electronic density, but they are primarily due to atomic 
stacking against the interface. The relative amplitude of 
the outermost oscillation increases with the valence and 
we attribute it to the drastic decrease undergone by the 
valence electronic density at the interface, which induces 
a steeper potential wall when moving from Li to Si. 



The self-consistent valence electronic density profiles 
are shown in figure H They oscillate near the surface 
although with a much smaller amplitude than the ionic 
ones. However, their relative phase exhibits an interest- 
ing behavior, which evolves from being in opposite phase 
for Li to almost in phase for Si. An opposite phase be- 
tween the ionic and the valence electronic oscillations had 
already been obtained in the Monte Carlo simulations of 
the interface of liquid alkalis^ and Ga^ (se e fig ures 7- 
9 in reference |j| and figure 11 in reference This 
behavior was attributed to the competition between the 
electronic kinetic energy contribution, which gets smaller 
values by weakening the oscillations, and the interaction 
term between electrons and ions, which being attractive 
takes smaller values for in-phase oscillations. An oppo- 
site phase was also found in OFAIMD simulations for 
the liquid-solid interface of Al^ and it was justified in 
terms of the interaction (represented by the use of a pseu- 
dopotential), between valence and core electrons, which 
tends to expel the valence electrons from the ionic posi- 
tions. Indeed, the idea that the ionic and valence elec- 
tron profiles oscillate in opposite phase appears widely 
accepted; however a closer scrutiny reveals some moot 
points: (i) the magnitude of the electronic kinetic energy 
is too small in comparison with the ion-electron interac- 
tion term (about 1%), so as to ascribe it a prime role in 
establishing the phase of the oscillations, and (ii) when 
(some of) the valence electrons are s-type a maximum of 
the valence electron pseudodensity is found at the ionic 
position (see the inset of figure EJ). Therefore, we have 
performed several tests to clarify the reasons underlying 
the phase-shift between the ionic and valence electronic 
density profiles. 

First a valence electron density profile was generated 
by superposing at ion sites pseudoatomic valence den- 
sities as obtained in the pseudopotential construction. 
This amounts to a linear response treatment of the va- 
lence electron density and therefore lacks any trace of 
the kind of competition argued above. These valence 
electron density profiles for Li and Si are compared with 
the self-consistent OFAIMD profiles in figure|3 and there 



is very good agreement. In particular the phase of the 
valence electron density oscillations is reproduced, which 
suggests that the phase difference between the ionic and 
valence electron density profiles is connected to some fea- 
ture of the pseudoatom density. The pseudoatom densi- 
ties projected onto the z-axis, are also shown in figure 
[3J The main features are the width and the presence of 
weak Friedel oscillations. To clarify the possible influ- 
ence of these features on the phase of the electron den- 
sity profile we have fitted the projected pseudoatom den- 
sity to a model with no Friedel oscillations. A good fit 
is obtained for a model density of the normalized form 
exp[— |z/cr| 3 ], which includes only a width parameter, a, 
with values: 1.60, 1.50, 1.29 and 1 .13 A for Li, Mg, Al 
and Si respectively. Superposing these model densities 
at the ionic positions generated by the simulations gives 
valence electron density profiles which, for all systems, 
are rather similar to the self-consistent OFAIMD ones, 
as shown in figure [3] for Li and Si. Again, the phase of 
the oscillations is preserved, and it is inferred that the 
Friedel oscillations in the valence pseudodensity are not 
responsible for the phase difference. The reason for the 
phase difference between ion and valence electron density 
profiles must lie in the width of the pseudoatom density 
as compared with the separation of layers in the ionic 
density profile. The ratio a/X has values 0.64, 0.58, 0.55 
and 0.45 for Li, Mg, Al and Si respectively, correlating 
with a decreasing phase difference between the ion and 
electron oscillations. Moreover, the OFAIMD results for 
other liquid metals near melting (Na, K, Rb, Cs, Ba, Ga, 
Tl, Sn) show a correlation between the phase difference 
and the ratio a/X. The systems fall into groups with sim- 
ilar phase differences: (i) the alkalis (0.62 < a/X < 0.64), 
(ii) Mg and Ba (0.58 and 0.59), (iii) Al (0.55) and (iv) 
Tl, Ga, Si and Sn (0.45 < a/X < 0.46). 

We stress that the width of the pseudoatomic density 
is characteristic of the atomic species, whereas the in- 
terlayer distance in the ionic profile, A, depends on the 
environment. For example, in the liquid-vapor interface 
of Al A = 2.35 A, whereas in contact with the (100) 
face of its solid fee phased Aw 2.1 A (which is close 
to the interlayer distance in the solid) and leads to a 
ratio a/X « 0.613, well within the range of the out-of- 
phase oscillations. To reinforce this argument we have 
taken the ionic positions of the liquid Al slab, and su- 
perposed the model pseudoatom densities with different 
widths, namely, a = 1.60, 1.29 and 1.00, with correspond- 
ing a/X = 0.68,0.55 and 0.43 respectively. Figure 
shows the resulting model valence electron density pro- 
files which evolve from opposite phase in the wider model 
to in-phase for the narrower one. 

The total electron density profile, which is the quantity 
accessible to experiment, is obtained by the sum of the 
self-consistent OFAIMD valence electron density profile 
plus the superposition of core electron densities, and is 
shown in figure|SJ Since the core densities are rather nar- 
row their superposition gives a profile in phase with the 
ion density profile. Consequently, when the valence elec- 
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FIG. 3: Valence electron density profiles at the liquid-vapor 
interfaces. The continuous line is the self-consistent result, 
the dashed line represents the linear superposition of pseu- 
doatomic densities, and the dotted line is the superposition 
of model densities without Friedel oscillations. The inset 
shows the projected pseudoatom electron densities (symbols) 
together with the fit to the model proposed in the text (lines). 
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FIG. 4: Normalized ion and valence electron densities of liquid 
Al obtained by superposition of several model densities of 
different widths, a. 



tron density is added the phase of the total electron den- 
sity profile depends on the relative weight of the core elec- 
tron (always in phase) and valence electron (any phase 
is possible) contributions. In this respect liquid Li is the 
most interesting case, since the valence contribution (in 
opposite phase) is 1/3 of the total while the core contri- 
bution (in phase) is 2/3. Figure shows that the core 
contribution dominates even for Li, and the total electron 
density profile is in phase with the ionic one. 

In summary, 20000 configurations of 2000-particle 
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FIG. 5: Normalized ion profile and total (core + valence) 
electron densities for the systems studied. The Mg, Al and Si 
data are shifted upwards by 0.5, 1 and 1.5 units respectively. 



slabs have been simulated ab initio using the OFAIMD 
method to obtain the ionic and valence electron density 
profiles of the liquid-vapor interfaces of Li, Mg, Al and 
Si; results have also been obtained and are reviewed here 
for some other systems. All the ionic profiles show layer- 
ing. The oscillations in the ionic profile are not induced 
by Friedel oscillations in the electron profile, but are due 
to atomic stacking. The valence electron density pro- 
files also show oscillations, but the phase with respect 
to the ion profiles evolves following a pattern that cor- 
relates directly with the ratio between the width of the 
pscudoatoms and the wavelength of the ionic oscillations. 
Nevertheless, the total electron density profile, even for 
Li, oscillates in phase with the ion profile, being domi- 
nated by the more localized and numerous core electrons. 
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